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ABSTRACT 

We present 3D implicit large eddy siraulations (ILES) of the turbulent convection in the envelope 
of a 5 Mq red giant star and in the oxygen-burning shell of a 23 Mq supernova progenitor. The 
numerical models are analyzed in the framework of ID Reynolds-Averaged Navier-Stokes (RANS) 
equations. The effects of pressure fluctuations are more important in the red giant model, owing to 
larger stratification of the convective zone. We show how this impacts different terms in the mean- 
field equations. We clarify the driving sources of kinetic energy, and show that the rate of turbulent 
dissipation is comparable to the convective luminosity. Although our flows have low Mach number 
and are nearly adiabatic, our analysis is general and can be applied to photospheric convection as 
well. The robustness of our analysis of turbulent convection is supported by the insensitivity of the 
mean-field balances to linear mesh resolution. We find robust results for the turbulent convection 
zone and the stable layers in the oxygen-burning shell model, and robust results everywhere in the 
red giant model, but the mean fields are not well converged in the narrow boundary regions (which 
contain steep gradients) in the oxygen-burning shell model. This last result illustrates the importance 
of unresolved physics at the convective boundary, which governs the mixing there. 



L INTRODUCTION 

Since the original publication of the rn ixing length 
theory (MLT) by iBohm-Vitensi (|T95l . much ef- 
fort ha s been devot ed to the improve ment of the 
theorv dCoughl [19771 FStcllingwcrf 19821 IXiond [1985 
i KuhfussI 119 86': Canuto fc Mazzitclh 1991; Canuto 1992"; 
Gehmevr fc Winklcr_ _1992; Wuchtcrl & Fcuclitinger 



1998tlDeng et al.li250(l . Nevertheless, MLT stih remains 
the standard choice in most state-of-the-art stellar 
evolution codes. 

With the wealth of data coming from asteroseismol- 
ogy missions (CoRot, Kepler), and expected from future 
observatories (Gaia, JWST), a new generation of stel- 
lar models is needed. The modeling of solar-like oscilla- 
tions requires reliable models for the Rey nolds stresses 
([Belkacem et all 120061: iSamadi et all I2012D . The inter- 
action between convection and pulsations, which sets 
the location of the red edge of the instability strip, 
needs a better time-dep endent theory of turbulent con- 
vection (jBuchler fc Koll ath 2000). In the deep interior, 
additional mixing is required at convective boundaries 
across the Hertzprung- Russell diagram, above convective 
cores (|Maedejll975t iMatraka et"al1ll982l : iSchroder et alj 
1997f). and belo w convective envelopes (Herwig 2000; 



Pace et aD l2012j ). Lacki ng a physical ly consistent de- 
scription of this process |Renzinilil987|) . extra-mixing is 
currently included in stellar evolution codes with ad-hoc 
parameterizations, so that predictive powers are ham- 
pered by the use of free parameters. Next generation 
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stellar evolution models should rely on a consistent de- 
scription of convective boundary mixing, together with 
the effect of internal waves induced by turbulence. 

The road to a satisfying theory of turbulent convection 
is difficult. Stellar convection is highly turbulent, with 
Reynolds and Rayleigh numbers having "astronomical" 
values (Re > 10^°, Ra > 10^°). The on-going develop- 
ment of computational physics allows numerical model- 
ing of turbulent systems having an increasing number 
of degrees of freedom, but at present no direct simula- 
tion of the problem is possible. Nevertheless, physical 
insight provided by computer simulations is invaluable 
in improving our understanding of stellar hydrodynam- 
ical processes. The path starting from large 3D data 
sets and ending with a recipe simple enough to be im- 
plemented in stellar evolution codes is not straightfor- 
ward. Ideally, a common framework should be used 
both for the analysis of multi-D data and for stellar 
evolution calculations, strengthening the underlying con- 
nection and making the projection from 3D to ID eas- 
ier. Reynolds- Averaged Navier-Stokes (RANS) equa- 
tions are a promising framework for this. Much effort 
has been already devoted to RANS in the co ntext of stel- 
lar hy d rodynamics, see e.g. [C anuto ( 199 71): iXiong et cd] 
([19971) :jCanuto fc Dubo^^dkovl (|1998tl : iDeng et all (|2006D : 
iCanutcl ()2011[ ) and references therein. We adopt the 
same methodology. 

A theory of turbulent convection should have a broad 
range of applicability: ideally to every type of star at all 
stages of evolution. Therefore, it is important to iden- 
tify fundamental properties of the physical process and 
to understand what changes with different stellar condi- 
tions. This is an ambitious project, which we start by 
considering two very different astrophysical cases: the 
convection in the envelope of a 5 Mq red giant star, and 
the convection in the oxygen- burning shell of a 23 Mq 
supernova progenitor. 

The paper is organized as follows: we present in Sect. [2] 
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Figure 1. Overview of the simulated region in the red giant model. Left panel: density (bottom curves) and temperature (top curves) 
stratification in the initial (ID) and 3D model rg.3D.mr. Right panel: squared Brunt- Vaisala frequency in the initial (ID) and 3D model 
rg.3D.mr. 

an overview of our stellar models and numerical tools. 
In Sect. |3l we develop a set of ID RANS equations as 
a framework to analyze our data. Results are presented 
in Sect. SI In Sect. [SJ we summarize our findings and 
conclude the paper. 



[20T1I) . To mimic surface cooling, we apply a Newtonian 
cooling term in the last 5 % of the star: 



T~Tn 



pCv 



fir), 



(1) 



2. STELLAR MODELS AND NUMERICAL METHODS 

In this section we describe the stellar models and the 
hydrodynamic codes, and summarize the properties of 
the 3D models. 

2.1. Red giant models 

The red giant model was presented in iViallet et al.l 
()2011|) . The initial ID structure was constructed by 
integrating the stellar structure equations with the fol- 
lowing input parameters: M = 5Mq, Tch = 4500 K, 
log(L/L0) = 3.1. In terms of stellar evolution phase, 
such a model would correspond to a 5 Mq star at the 
end of central He burning, finishing its blue loop and 
evolving toward lower T^ff and awa y from the red edge 
of the Cepheid Instability Strip (see lAlibert et al1ll999l 
their Fig. 1 and Tab. 6). The stellar structure code is de- 
scribed in Bara fFe fc El Eid (19911 . It uses mixing- length 
theory (with a = 1.7) to treat convection, and the extent 
of the convective region is based on the Schwarzschild 
criterion. The structure was integrated from the photo- 
sphere down to 20 % of the stellar radius, stopping to 
avoid the nuclear burning region. This initial stratifica- 
tion is used as an input model for the multi-D hydro- 
dynamic code. The initial stellar structure is shown in 
Fig. [TJ The left panel shows the temperature and den- 
sity stratifications. The model is characterized by a total 
density stratification log(/9bottom/ptop) 4.4. The total 
pressure stratification (not shown) is log(pbottom/-Ptop) ~ 
6.2, or equivalently ^ 14.3 pressure scale-heights. The 
right panel shows the radial profile of N'^, the Brunt- 
Vaisala frequency squared. The convective region ex- 
tends down to r ^ 2.3 x IQ-'^^ cm, nearly half of the 
star in radius. The surface layers are characterized by a 
strong s uperadiabatic strat ification. 

As in IViallet etall (poTTI ). we use a proxy for the sur- 
face layers. The surface layers are numerically difficult to 
handle: the decrease in the pressure scale-height yields 
small scale convective eddies which are difficult to resolve 
when a single grid is used (see discussion in iViallet et al.i 



where /(r) is a spatially varying function that is equal 
to 1 above a given radius Vc with a smooth transition to 
zero, T is the cooling timescale and Tq is the forci ng tem- 
peratu re. We use the same parameters as iu iViallet et al.l 



(|20TTI) : 



0.95i?^, To = 32 750 K, and r 



10^ 



About 4.5 pressure scale-heights of the initial ID model 
are absorbed into the Newtonian cooling region. 

We solve the equations describing the evolution of den- 
sity, momentum, and total energy for a single fluid, tak- 
ing into account gravity, radiative diffusion, and surface 
cooling: 



d 



---V-{pu), (2) 
: -V • {pit ®u)-Vp + pg, (3) 
: -V • {petu + pu)+ pu-g + V ■ (xVT) - ^,(4) 



where p is the density, u the velocity, et = f-i + f-k the 
specific total energy (e^ is the specific internal energy 
and Efe the specific kinetic energy), p the gas pressure, 
T the temperature, g the gravitational acceleration, and 
X the thermal conductivity. For photons, the thermal 
conductivity is given by 



X 



3kp ' 



(5) 



where k is the Rosseland mean opacity, and a the Stcfan- 
Boltzmann constant. 

The 3D simulations were performed w it h the MU- 
SIC code, a s described in IViallet et al.l (j201lD and 
IViallet et ahl (|2013) (submitted to A&A). The spatial 
method is based on a finite volume discretization of a 
spherical wedge ( "box in a star" approach) . The method 
is based on staggered velocity components, and has been 
extended to 3D for this work. The time-marching scheme 
used for the models presented in this paper is based 
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Figure 2. Left panel: evolution of the total kinetic energy in the i 
energy in model rg.3D.mr. 

on the Mi nimum Re s idual Approximate Implicit Scheme 
from Bot chev et al.l ()1999[) . The MUSIC code is opti- 
mized to run on parallel computers and it uses domain- 
decomposition to distribute the computation over the 
computational nodes. The Message Passing Interface 
(MPI) is used to handle communications of boundary 
data. We use periodic boundary conditions in horizontal 
directions, and non-penetrative stress-free conditions at 
the bottom and top of the domain. As the nuclear burn- 
ing region is not included in the computational domain, 
a radiative flux corresponding to the stellar luminosity is 
imposed at the inner boundary. 

We do not model viscosity explicitly in the equations. 
The expected value of the molecular viscosity in stel- 
lar interiors implies huge Reynolds numbers (larger than 
10^°). It is therefore impossible to model all scales of 
the flow, from the stellar scale down to the dissipation 
scale, on current generation of computers. We adopt 
the Implicit Large E ddy Simulation paradigm (ILES 
iGrinstein et al.l I2007D and solve the inviscid equations 
to model the turbulent flow. The underlying motivation 
of ILES is that monotonic, finite- volume based methods 
have physical properties of the Navier-Stokes equations 
"built-in" within the numerics (unlike spectral and fi- 
nite difference methods). The conservation properties 
of finite volumes schemes and the monotonicity preserv- 
ing property enforce correct physical behavior at the grid 
scale. As a result, the loss of information that takes place 
at the grid scale mimics turbulent dissipation (see l4.1.2|) . 

Models with two different resolutions were computed; 
model rg.3D.lr with a 216 x 128^ resolution, and model 
rg.3D.mr with a 432 x 256^ resolution. The main proper- 
ties of these models are summarized in Table [T] We use 
an opening angle of 45° X 45° . Model rg. 3D. I r was evolved 
for 3100 days of model time (for ^ 3.6 x 10"* CPU hours on 
512 cores) and model rg.3D.mr was evolved for 1050 days 
of model time (for ~ 2 x 10^ CPU hours on 2048 cores), 
starting from the solution a.tt = 2000 d of model rg.3D.lr. 
The computations were performed on the IBM iDataPlex 
supercomputer "Hydra" at the Rechen-zentrum Garch- 
ing, and on the SGI Altix ICE 8200 supercomputer "Zen" 
at the University of Exeter. At the beginning of the com- 
putation, the Newtonian cooling term modifies the top 
layers and this triggers the convective instability. Down- 
drafts form and sink. They break up rapidly due to hy- 
drodynamical instabilities, and the fiow becomes turbu- 
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giant models. Right panel: space-time diagram of the specific kinetic 
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Figure 3. Radial profiles of the rms fluctuations of density, tem- 
perature, and pressure (continuous lines) in model rg.3D.mr. The 
black das hed line shows pu'"^ /p for comparison (see discussion in 
Sect. [431 . 

lent. Left panel of Fig. [2]shows the evolution of the total 
kinetic energy in the domain. The initial rise of the total 
kinetic energy is strong, followed by a slow decay toward 
a quasi steady state. The left panel of Fig. [1] shows the 
temperature and density stratifications obtained in the 
3D model when a quasi-steady state is reached. The ef- 
fect of the Newtonian cooling is to drive an isothermal 
region at the top. The right panel of Fig. [T] shows the 
profile of N"^ in the 3D model. The most striking fea- 
ture is the modification of the radial profile of N"^ , which 
shows a much steeper slope at the bottom boundary of 
the convective zone in the 3D model than in the initial ID 
model. This is an evidence for overshooting, and will be 
discussed in Sect. 14.61 In the 3D model, the density strat- 
ification is log(/5bottoin/ptop) ^ 2.3 in the convective zone. 
The pressure stratification is log(pbottom/ptop) ~ 3.4 in 
the convective zone, or equivalently ~ 7.8 pressure scale- 
heights. Clarifying the influence of such a large stratifi- 
cation on the properties of the turbulent convection is an 
important focus of this work. Figure [3] shows the radial 
profiles of the rms values of the density, temperature, and 
pressure fiuctuations, normalized by their mean values. 
In the bulk of the convective region, all fluctuations have 
the same relative order of magnitude. 

Computer animations and snapshots of the models 
illustrate the strong asymmetry of the flow, with the 
presence of plumes triggered by cooling at the sur- 
face that sink in the convective zone within much 
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Figure 4. Left panel: snapshots of the computational domain for model rg.3d.mr. The magnitude of the vorticity is shown (two 
perpendicular cuts in the vertical direction were taken and are plotted together). Right panel: Radial profiles of rms velocity components 
for model rg.3d.mr : Ur (thick solid line), ug (thin solid line), and (thin dashed line). 



slower upfLowing material. This is 
gree of strat ification, as empha sized 



due to the de- 
by previous SD 

simulations (IStein fc N ordlund Il989t i Cattaneo et al 
19911: IBrummell et all 11996: .Porter fc Woodward .200 " 



Brummell et alJ 120021 ). These plumes form a "network" 



of downflows in strong interaction, and from time to time 
plumes coalesce to form a strong downdraft which sinks 
through the whole stratification. The right panel of Fig. 
[2] shows a space-time diagram of the specific kinetic en- 
ergy, which shows the imprint of these sporadic events. 
The left panel of Fig. |4] is a snapshot of the flow, as seen 
in the magnitude of the vorticity vector field uj = V x u 
and emphasizes the prominence of small scale structures, 
characteristic of 3D turbulence. The right panel of Fig. 
[3] shows the rms values of the velocity components Ur, 
Ug, u^. The tangential and azimuthal components are 
roughly equal as the angular directions are homogeneous. 
Including rotation and/or magnetic fields would break 
this symmetry. The corresponding Mach number goes 
from ~ 0.1 at the top of the convective zo ne to ~ 0.01 
at the bottom. As in lArnett et al.l (|2009r ). we define a 
global rms velocity such as 



1 



2 

rms 



k,CZ, 



(6) 



where Mcz and -Ek.cz are the total mass and kinetic 
energy in the convective zone (CZ), respectively. The 
convective zone is taken as the region between rin = 2 x 
10^^ cm and rin = 4 x 10^^ cm, but we have checked that 
results are not too sensitive to these values. We find 
I'rms = 2.34 X 10^ cm/s for model rg.3D.mr. We define 
the convective turn-over timescale as 

r, ^cz 



(7) 



which yields Tconv = 198 d for model rg.3D.mr. 

Finally, Fig. 0] and the right panel of Fig. [5] hint at 
the presence of g-modes in the underlying radiative zone. 
The modes are excited at the convective boundary layer. 
In this paper, we will focus on the dynamics in the con- 



Table 1 

Summary of the Red Giant Simulations. 



Parameter 


rg.3D.lr 


rg.3D.mr 


Grid zoning 


216x128^ 


432x256^ 


nn, ''out (10^^ cm) 


0.82, 4.09 


0.82, 4.09 


A6», A0 


45°, 45° 


45°, 45° 


CZ stratification {Hp) 


7.8 


7.8 


iav(Atav) (days) 


2650(800) 


2650(800) 


frms (10^ Cm/s) 


2.27 


2.34 


Tconv (days) 


204 


198 


(10^'^ erg/s) 


4.9 


4.9 


Ld (10^6 erg/s) 


7.33 


7.39 


/d (10" cm) 


7.0 


7.7 


Td (days) 


18 


19 


Pe 


4900 


5200 



Note. 

^rms : global rms velocity, Tconv : 

convective turnover timescale, L*: luminosity of 
the stellar model, L^: rate of kinetic energy dis- 
sipation, l^: dissipation length-scale, t^: dissipa- 
tion time-scale, Pe: Peclet number. 

vective region, and leave the study of wave excitation 
and dynamics for future work. 

2.2. Oxygen-burning shell models 

We simulate the turbulent flow in a convective oxygen 
burning shell using as initial conditions a 23 Mq stellar 
model of solar composition that was previ ously evolved 
with the TYCHO stellar evolution code (|Arnett et al.l 
|2009() up to a point where oxygen, neon, carbon, helium, 
and hydrogen are burning in concentric shells about a 
silicon-sulfur rich c ore. Additio nal det ails can be found in 
iMeakin fc ArnettI ()2006l I2007D . As in iMeakin fc ArnettI 
( 20071) . we study only the oxygen burning convection 
zone and its stably stratifled bounding layers. 

Reactive-hydrodyna mic evolution was simula ted with 
the PROMPI code (jMeakin fc ArnettI [2007I) . a dis- 
tribut ed memory adapt ation of the PROMETHEUS 
code (|Frvxell et al.l [19891) . an Eulerian implementation 
of the piec ewise parabolic method (PP M) hydrodynam- 
ics scheme (jColella fc Woodw ard"1984') extended to treat 
realistic equations of state (jColella fc Glaz.l985i ). multi- 
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Figure 5. Overview of the simulated region in the oxygen-burning shell model for the initial model (ID) and for the 3D model ob.3D.mr. 
Left panel: density (continuous line) and temperature (dashed line) stratification. Right panel: squared Brunt- Vaisala frequency. 



species advection, and a general nuclear reaction net- 
work. The setup for th ese simulations are nearl y identi- 
cal to that described in lMeakin fc Arnettl ()20Q7[ ). includ- 
ing the 25 species followed to track nuclear evolution, 
and differ only in terms of resolution and domain size. 
As with the MUSIC code, PROMPI solves the inviscid 
equations, and we rely on the ILES paradigm to model 
sub-grid scale dissipation. The major differences between 
the two is that PROMPI is multi-fluid, includes nuclear 
burning, and that rad iative diffusion is negligible in this 
context ([Arnettl [l996f l. Figure [5] shows the initial den- 
sity and temperature stratifications (left panel), and the 
initial profile of A'^^ (right panel). The latter is charac- 
terized by a narrow peak at the bottom boundary of the 
convective zone, due to the sharp composition gradient 
that characterizes the boundar y of the nuclear burning 
region. We refer the reader to (Meakin fc Arnet'^ ()2007[ ) 
for similar figures as those shown earlier for the red giant 
model. 

As with the red giant models, the oxygen shell burn- 
ing models studied in this paper use a 45° x 45° open- 
ing angle. Models with three different spatial resolu- 
tions were computed: model ob.3D.lr with 192 x 128^ 
zones, model ob.3D.mr with 384 x 256^ zones, and model 
ob.3D.hr with 768 x 512^ zones. The main properties of 
these models are summarized in Table [2j Models ob.3D.lr 
and ob.3D.mr are used for comparison with the red giant 
models, whereas model ob.3D.hr is used to further assess 
the numerical convergence of our results with resolution 
(see Sect. 14. 7|) . All of the oxygen burning models were 
computed on the University of Tennessee's Kraken Cray 
XT5. Model ob.3D.lr was evolved for approximately 600 s 
of model time (on 768 cores for ~ 5 x IC* CPU-hours). 
Model ob.3D.mr was evolved for approximately 280 s of 
model time (on 12,288 cores for - 3.5 x 10^ CPU-hours) 
starting from the solution at i = 300 s of model ob.3D.lr. 
Model ob.3D.hr was evolved for approximately 200 s of 
model time (on 24,576 cores for 4 x 10® CPU-hours) 
starting from the solution a.tt = 310 s of model ob.3D.mr. 

At the beginning of the computation, the convective 
instability is triggered by a band of random, small am- 
plitude density perturbations {p' / p ^ lO^'') imposed 
within the convection zone. The flow rapidly becomes 
turbulent as it fills the convectively unstable region. The 
early transient evolution of the models is characterized 
by a strong penetration of the fiow in the top stable re- 
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Figure 6. Same as Fig. [S] but for the oxygen-burning shell model 
ob.3D.mr. 

gion. After roughly 300 s, a quasi steady-state obtains 
with a slow evolution resulting from a net heating due 
to the imbalance between nuclear burning and neutrino 
cooling (a common feature of neutrino-cooled stages of 
stellar evolution which results in growth of the convec- 
tive region). As a result, global characteristics of the 
model (e.g., total nuclear burning, total kinetic energy, 
etc.) increase slowly with time. The nuclear evolution 
time scale for this phase is roughly 5 x 10"^ s. The 
convective region is characterized by a density stratifi- 
cation Pbottom/ptop ~ 6, and a pressure stratification 
Pbottom/ptop ~ 7.5 or two pressure scale-heights. The 
global rms velocity and turnover timescale are computed 
as for the red giant (see previous section), and shown in 
Table [2] Figure |6] shows the radial profiles of the rms 
values of the density, temperature, and pressure fluctua- 
tions, normalized by their averaged values. The largest 
values in the relative magnitude of the fluctuations occur 
at the boundaries, and are of the order of a percent. 

3. HORIZONTALLY-AVERAGED MEAN-FIELD EQUATIONS 

Using the RANS methodology, we develop in this sec- 
tion a framework of mean-field equations which are con- 
sistent with spherical geometry. A similar approach 
has been appl ie d to the k in etic energy equat i on; se e 
IHurlburt et ajl (IT9861 IT991 : IMeakin fc A^^ ([20?)1 . 
We show here how this can be extended, as the ki- 
netic energy is only one of the several mean-field equa- 
tions that can be derived from the hydrodynamical equa- 
tions. Although we refer to our equations as "Reynolds- 
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Table 2 

Summary of the Oxygen Burning Simulations. 



Parameter 


ob.3D.lr 


ob.3D . mr 


ob.3D.hr 


Grid zoning 


192x128^ 


384x256^ 


786x512^ 


rin, rout (lO'-* cm) 


0.3, 1.0 


0.3, 1.0 


0.3, 1.0 


Ae, A(f> 


45°, 45° 


45°, 45° 


45°, 45° 


CZ stratification {Hp) 


2 


2 


2 


tav{Atav) (s) 


494(230) 


494(230) 


429(165) 


■"rms (10® cm/s) 


9.2 


9.57 


9.2 


'^conv (s) 


96 


92 


96 


Lnuc (10*6 erg/s) 


2.69 


2.63 


2.47 


Ld (10^6 erg/s) 


0.30 


0.29 


0.30 


«d (10* cm) 


4.9 


5.5 


4.7 


Td (s) 


26 


29 


26 



Note. — frms; global rms velocity, rconv: convective 
turnover timescale, Lnuc^ total energy released by nuclear 
burning, L^: rate of kinetic energy dissipation, i^; dissipa- 
tion length-scale, t^: dissipation time-scale. 



Averaged Navier-Stokes" equations, we actually intro- 
duce the effect of viscosity only through the kinetic en- 
ergy dissipation rate (units: erg s""'^ g""'^)- In the 
large Reynolds number regime, the viscous dissipation is 
the dominant contribution of viscosity in the mean-field 
equations. This is the so-called "dissipation anomaly" 
characteristic of 3D turbulence, see discussion in Sect. 
14.71 Other terms, such as the viscous flux of kinetic en- 
ergy, are neglected. 

To obtain our set of ID equations, we introduce two 
types of averaging: 1) statistical averaging, denoted by 
an overbar; 2) horizontal averaging, denoted by brack- 
ets. Practically, statistical averages are computed by 
performing a time average (the ergodic hypothesis). The 
horizontal average of quantity q is defined as: 

{q) = ^J^j{r,e,^)dn, (8) 

where dQ — sinOdOdcf) is the solid angle in spherical co- 
ordinates. 

We decompose the fiow quantities into mean and fluc- 
tuating components: 

q = ¥)+q, (9) 

so that (g') = 0, by construction. We also introduce 
Favre averaging, which is a density-weighted average: 



- % (10) 
{p) 

and which leads to the decomposition 

q={q)+q"- (H) 

The density and energy equations lead to their av- 
eraged counter-part. From the momentum equations, 
the three mean-fleld equations which are consistent with 

spherical geometry concern the mean radial velocity {ur), 
the mean specific angular momentum along the z-axis 

(jz), and the mean specific kinetic energy {tk)- By "con- 
sistent with spherical geometry" we mean that the re- 
sulting equations do not have terms which have an ex- 
plicit angular dependance, as would for instance have the 



mean-field equations for ug or u^. The energy equation 
can be formulated either for the mean specific internal 

energy (e^), the mean specific total energy (ej), or the 

mean specific entropy (s). For completeness, we show 
all forms. The resulting equations are summarized in 

Table [3] We show the equation for {jz) although it is 
"trivial" when rotation is not included; we will not dis- 
cuss it further here. Note, however, that when rotation 
is included, the horizontal average as introduced here 
is not suited because the angular directions are not ho- 
mogeneous anymore. In this case, similar ID averaged 
equations cannot be obtained. 

We denote by the radial part of the divergence op- 
erator in spherical coordinates, i.e. V^/ = -pidr{r^f)- 
The mean-field equations are characterized by second- 
order correlations, which stem from the Reynolds/Favre 
decompositions, as, e.g., the turbulent fluxes which are of 

the form {q"u'^) or (q'u'j.). The equations in Table |3] are 
written in terms of the averaged Lagrangian derivative 

{D^)q = dtq + Mdrq. (12) 

The connection with the Eulerian, conservative form is 
immediate: 

Jp}{D^)q = dt({p}q) + Vr{qM), (13) 

where we used the ID averaged continuity equation. 

Connection with the common form of the stellar evo- 
lution equations can be made using: 

dr = ^Trr^Jp)dm, (28) 
{Ih)^dt\m (29) 

The flrst equation is the standard relation between 
the Eulerian derivative at constant radius and the La- 
grangian derivative at constant mass shell. The sec- 
ond equation states that (Dt) is exactly the Lagrangian 
derivative at constant mass. The resulting equations are 
summarized in Table S) 

4. RESULTS 

Unless stated otherwise, the results in this section are 
based on models rg.3D.mr and ob.3D.mr. In Sect. 14. ll we 
use the mean-field equations developed in Sect.[3]to ana- 
lyze our 3D data. These ID equations provide a powerful 
framework to analyze the physical processes characteriz- 
ing the numerical models, and to assess the consistency 
of the models with the physical equations. Section I4.2l is 
devoted to the analysis of the turbulent velocity field. We 
motivate the distinction between "deep convection" , rel- 
evant for the red giant models, and "shallow convection" , 
relevant for the oxygen-burning shell models. We analyze 
the anisotropy of the Reynolds stresses, which charac- 
terizes the asymmetry of the flow, and which motivates 
a decomposition of the velocity fleld into a large scale 
component, characterizing the plumes, and an isotropic 
component at small scales. We show that the isotropic 
component provides a natural interpretation of the ki- 
netic energy damping in terms of a dissipation in a tur- 
bulent cascade. In Sect. 14.31 we show how the magnitude 
of pressure fluctuations is related to the stratification of 
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Table 3 

ID RANS equations in Lagrangian form. 



V)(Pt){^) 
Tp){Dt)l^) 

W){Dt)(^) 



(p)Vr(Mr> 



■ SJr(p)(t'{K) + Vr{x)dr{T) - {p)V r (Ur) + Vr{x'drT') - {p'V ■ u') + {ptnuc) + (pea) 



■ V,(p>(/i"0 - V.(p>(e''0 - (p>V,K> + V,(x>ar(T> + Wrix'drT') + (penuc) 



Vrjp}{^;) - (^V . Fr) + ip'-^^^^) 



V,(p>(e''<> - V.(p'mO + (P'V • «') + {p'u' ■ g) ~ (pea) 



: - V,.(p>(<^ - dr{p) - {p){g) + ^(2(efe> - K)' 

r 



(14) 

(15) 
(16) 

(17) 

(18) 

(19) 

(20) 



Note. — Definitions: density p, temperature T, pressure p, velocity u, radial velocity component Ur, specific internal energy ti, specific 
kinetic energy ej., specific total energy et, specific entropy s, 2-component of the specific angular momentum = rsindu^, thermal 

conductivity Xi radiative flux Fr, gravitational acceleration g, rate of nuclear energy production Enuc, rate of viscous dissipation e^. 

Reynolds decomposition: q = {q) + q' , Favre decomposition: q = {q) + q" ■ 

Table 4 

ID RANS equations in Lagrangian mass coordinate (same deflnitions as in Table [HJ. 



dtr 

dt{Ur) 

dt{j7) 



(21) 



- 9„(47rr2(p)(^)) + dm{'inr^{x)dr{T)) ~ (p>a™(47rr2(«,)) + dmi^Trr^ {x'drT')) - + {^) + (e^) (22) 

- e„(47rr2"(^(h^>) - dm{47Tr^Jp)(^)) + dmiinr^'i^drlX)) + dm(4TTr^ {x' drT')) - l^dmiiTTr^i^)) + (^) (23) 

- 9™(4^r2(^(5x'» + ■ Fr) + i^(p ^""'- + ^'i ) (24) 



- d^{4nr^{p)(^)) ~ e„(4^r2(p'<» + + ^^^^ - (e,) 

(P> (p) 

1 



„ =-a„(4vrr2^(<^» 



the convcctive zone, which explains the main differences 
in the mean field analysis between the red giant model 
and the oxygen-burning shell model. In Sect. WA[ we an- 
alyze the turbulent fluxes. These are the quantities that 
should be modeled in a theory of turbulent convection. 
In Sect. 14.51 we discuss the kinetic energy driving and the 
turbulent dissipation. We show that in steady state the 
rate of turbulent dissipation is of the order of the convcc- 
tive luminosity. Section 14.61 is devoted to the red giant 
models, for which we discuss thermal effects and charac- 
terize the overshooting at the bottom of the convective 
region. Finally, we discuss in Sect. I4.7l the convergence of 
our results with resolution, and consider the implications 
regarding the global dynamics in the numerical models. 

4.1. Mean-field analysis 

We use the framework developed in Sect. [3] to ana- 
lyze the data from our multi-D simulations, in terms of 
the radial momentum, kinetic energy, and total energy 
balances. The temporal averaging is performed on the 
interval [tav — Aiav/2,iav + Atav/2], with tav and Aiav 
given in Tables [T] and for each model. In practice, we 
find that averaging the data over more than two turn- 
over timescales yields very similar results. In this section. 



(25) 

(26) 
(27) 



Table 5 

Radial momentum balance for model rg.3D.mr - integral 
budget over the convective zone (g.cm.s~^). 



Term 


Value 


Term 


Value 




3.85(28) 




5.62(30) 


-/V.(p)(<^dV 


6.03(28) 


Residual 


-5.91(28) 


-n9r{p)-{p){g))dV 


-5.66(30) 







the models were averaged over the common time period 
available at different resolutions, namely four turnover 
timescales (800 d) for the red giant model, and roughly 
2.5 turnover timescales (230 s) for the oxygen-burning 
shell model. 

Throughout the paper, we use the lower case letter / 

to denote turbulent fluxes, / — {p){urq") or f — {u'^q'), 
depending on the quantity q. 

4.1.1. Radial momentum balance 
Equation ([19)) can be written as: 
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Figure 7. Mean-field analysis for model rg.3D.mr, time averaging is performed over 800 days. From top to bottom: radial expansion 
velocity, kinetic energy, and total energy balance. Left panel: Radial profiles of the terms in the balance. For the sake of clarity, each 
term is multiplied by 47rr^. As a result, the volume integral of a term is equal to the area below the corresponding curve, and it can be 
appreciated visually. Right panel: Integral budget over the convective zone. 



Table 6 

Radial momentum balance for model ob.3D.mr - integral 
budget over the convective zone (g.cm.s~^). 



Term Value Term Value 



-fJ^Dti^V -4.55(31) JJ^:i!^dV 1.88(38) 
- / Vr- (p>"{<^V -2.13(36) Residual -4.39(37) 
- J (drW) -W){9))<iV -1.42(38) 



{p){dt){Ur) = - Wr{p)K') - dr{p) - (p) (<?) 

+ (30) 
r 

This equation is an equation for the mean radial compo- 
nent of the mass flux since by definition (p){ur) — (pur). 
The right-hand side involves the divergence of the ra- 
dial component of the Reynolds stress Rm the bal- 
ance between the gradient of the mean pressure and the 
mean gravity force, and a geometric term characteristic 
of spherical geometry in which a horizontal velocity Vh 
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Figure 8. Same as Fig. [7]for model ob.3D.mr, with time averaging performed over 230 s. 



induces a radial acceleration v\lr due to inertia. 

The top-left panels in Fig. [7]and Fig. |8]sliow the radial 
profiles of these terms in models rg.3D.mr and ob.3D.mr. 
In both cases, we find the Lagrangian time derivative to 
be negligible, showing that the balance is in a statistically 
steady state. The figures show that a slight hydrostatic 
imbalance is due to the turbulent ram-pressure, with a 
smaller, but not negligible, contribution from the inertial 
acceleration due to horizontal motions. Defining 



Pturb = 



rb/p decreases from 

-4 



(31) 
2 X 10^2 



we find that the ratio ptu 

at the top of the CZ to ^ 2 x 10^"' at the bottom m 
the red giant model. In the oxygen burning-shell model, 
Pturb/p ~ 10"'^ in the bulk of the convective zone, and it 



decreases rapidly at the boundaries. 

In the red giant model, the hydrostatic imbalance cor- 
responds to an inward acceleration of ~ 3.5 % g at the 
top of the convective zone, and an upward acceleration of 
^ 0.05 % at the bottom. For the oxygen-burning shell, 
the results are similar, with values within ±0.2 % oi g 
throughout the convective zone. Thus, as expected for 
the deep interior, hydrostatic equilibrium holds, typically 
to the order of a percent or less. In the photospheric re- 
gions the turbulent pressure can significantl y affect the 
hydrostatic balance (jStein fc Nordlundlll989D . 

The top-right panels in Fig. [7] and Fig. |S] (see also 
Tables [5] and [6]) show the balance upon integration over 
the convective zone. The residual is very small in the red 
giant models, showing good consistency with the phys- 
ical equations. The oxygen-burning shell models show 
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Table 7 

Kinetic energy balance for model rg.3D.mr - integral 
budget over the convective zone (erg s~^). 



Term 


Value 


Term 


Value 


-JVrfkdV 
- f \7rfpdV 


-3.47(33) 
-1.61(32) 
-2.68(34) - 


JWpdV 
JWtdV 
f{P^d)dV 


2.48(36) 
4.94(36) 
-7.39(36) 


Table 8 

Kinetic energy balance for model ob.3D.mr - 
budget over the convective zone (erg s" 


integral 


Term 


Value 


Term 


Value 


-/I^Dt(lI>dy 
-JVrfkdV 

- f ^rfpdV 


-1.91(44) 
-1.96(41) 
-1.04(42) - 


JWpdV 
JWidV 
-f{P^d)dV 


-1.45(44) 
3.27(45) 
-2.93(45) 



spurious oscillations at the inner boundary, whicli are 
responsible for the non-zero residual seen in the figure. 
This is due to the steep gradients at this boundary (see 
Sect. I4.7p . and we have checked that the consistency is 
very good everywhere else. Prog ress has been made in 
dealing with multi-fluids in PPM (jPlewa fc Mulledll999l : 
IWoodward et al. 2008), which may help. 

4.1.2. Kinetic energy balance 

At the first sight, the temporal behavior of the kinetic 
energy (KE) shows a high level of complexity (see Fig- 
ure [2] for the red giant model). The mean kinetic en- 
ergy balance, Eq. ([TSl) . allo ws insight in the dynamics 
of the turbulent convec tion (jHurlburt et al.lll986l [l99l 
iMeakin fc ArnettI |2007[) . The KE equation involves two 

turbulent fluxes: the KE turbulent flux fk = {p){u'^e'l) 

and the acoustic flux fp = (p'u'j,). These transport 
terms characterize the non-locality of turbulent convec- 
tion. Source terms for KE are the gravitational work 
due to density fluctuations Wb — {p'u' ■ g) , and the work 

done by pressure fluctuations Wp — (p'V • m'), also called 
"pressure-dilatation" . There is some freedom in the for- 
mulation of the kinetic energy equation, especially in the 
expression /splitting of the drivin g te rms, see e.g. dis- 
cussions m iNordlund et~all (|2009| ) and lArnett fc MeakinI 
(pOlii) . Here we choose to split the driving into Wp and 
Wh as these are thermodynamically the most fundamen- 
tal quantities: Wb measures the reversible exchange be- 
tween kinetic energy and potential energy, Wp measures 
the reversible exchange between kinetic energy and in- 
ternal energy. Furthermore, this splitting will highlight 
the differences between the red giant and oxygen-burning 
shell models. Wb is often called the buoyancy work, an 
appellation inherited from the Boussinesq literature. In 
Sect. 14.51 we present another expression which empha- 
sizes the underlying physical origin of the driving. 

As discussed in Sect. [21 in the ILES paradigm we rely 
on the numerical dissipation at the grid scale to mimic 
the effect of viscosity. As a result, kinetic energy is dis- 
sipated at the grid scale, and a sink term —ptd is in- 
troduced in the analysis. We compute each term in the 
balance equation, and identify the residual with the dis- 
sipation. This is valid for conservative methods of ad- 



Table 9 

Total energy balance for model rg.3D.mr - integral budget 
over the convective zone (erg s^^). 



Term 


Value 


Term 


Value 


-!V)Dti^)dV 


2.18(36) 


-jW}\7r{^)dV 


-2.13(36) 


^JVrhdV 


-1.01(35) 


-JqdV 


-5.45(36) 


-JVrfkdV 


-1.61(32) 


Residual 


6.30(35) 


-f\'r(Fr)dV 


4.88(36) 







vection. 

The middle-left panels of Fig. [7] and Fig. [S] show 
the radial profiles of the different terms in the KE bal- 
ance for models rg.3D.mr and ob.3D.mr. The inferred 
dissipation is shown as a black dashed line. The ki- 
netic energy balance is in a statistically steady state, as 
the contribution from the time derivative is negligible. 
In the red giant model, we find that both Wb and Wp 
drive (i.e. positive contribution) turbulent motion within 
the convective zone. Wb dominates the driving but the 
pressure-dilatation term is significant and cannot be ne- 
glected. This contrasts with the oxygen-burning shell 
models where driving is largely dominated by Wb, with 
no significant contribution from pressure-dilatation. We 
attribute the importance of Wp in the red giant model 
to the degree of stratification, which is larger than in 
oxygen-burning model (see Sects. 14.31 and [4.5p . The con- 
vective region is bounded by regions where Wb < 0. The 
KE dissipation profile has a large amplitude and is rather 
evenly distributed over the convective region. There is 
no local balance between the driving and the dissipa- 
tion, due to transport by the turbulent fluxes. In the 
red giant model, the net effect of the turbulent KE flux 
is to transport energy downward in the convective zone, 
whereas in the oxygen-burning shell model KE is trans- 
ported upward. In the red giant model, the acoustic flux 
has an opposite effect to the KE flux as it transports 
KE upwards. It has a more complex behavior close to 
the boundaries, due to an important contribution from 
waves. In the oxygen-burning shell model the acoustic 
flux is important close to the boundaries, but it has a 
negligible effect in the bulk of the convective zone. 

The middle-right panel of Fig. [7] and Fig. [5] (see also 
Tables [7] and [8]) show the KE integral budget over the 
convective zone. Upon integration on the convective 
zone, we are left with a balance between driving and 
dissipation: 

J WbdV + J WpdV « J J^dV. (32) 

To obtain the balance described by Eq. (|32l) . we have 
integrated over the unstable region and the convective 
boundary layers, where dissipation is taking place. The 
rough balance between driving and dissipation obtained 
over this region implies that only a small amount of ki- 
netic energy in transmitted to the stable zones in form 
of internal waves, so that the dissipation of waves energy 
in these regions is small. 

4.1.3. Total energy balance 

We now discuss the total energy balance, Eq. (|T6)) . 
The fluxes that appear in this equation are the radiative 
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Figure 9. Radial profiles of different luminosities in rg.3D.mr (left panel) and in ob.3D.mr (right panel). 



Table 10 

Total energy balance for model ob.3d.mr - integral budget 
over the convcctivc zone (erg s^^). 



Term 


Value 


Term 


Value 


-jW)Dt{^)dV 


-1.92(46) 


-/"(pyvrK)dv 


-8.20(45) 


-/V./hdV 


4.26(43) 


- / (pi;nuc>dV 


2.73(46) 


-JVrfkdV 


-1.96(41) 


Residual 


1.20(44) 



flux, the enthalpy flux, and the kinetic energy flux: 



(33) 
(34) 
(35) 



We discuss this form of the energy equation because 
of its relevance for stellar evolution calculations. For in- 
stance, let us consider the simplest case by assuming a 

steady state and (ur) = 0. Multiplying Eq. ([TB]) by Anr^ 
and integrating over radius, one obtains: 



(36) 



where L{r) — {p€nuc)dV. Stellar evolution codes use 
the mixing-length theory to compute fh and ignore fk- 
Figure [9] illustrates the different terms in Eq. ([36|) for 
models rg.3D.mr and ob.3D.mr. The radiative flux is neg- 
ligible in the oxygen-burning shell models and is there- 
fore ignored. In the red giant model, the radiative lu- 
minosity is equal to the stellar luminosity in the radia- 
tive zone, and it decreases in the convective zone where 
a large and outward directed enthalpy flux takes over. 
Furthermore, the red giant model is characterized by a 
downward directed kinetic energy flux reaching an am- 
plitude of roughly 35 % of the maximum enthalpy flux. 
The convective boundary is characterized by a negative, 
i.e. downward directed, enthalpy flux. This is somewhat 
counterbalanced by a bump in the radiative luminosity at 
the same location. The total luminosity is not constant, 
which means that the model is not in thermal equilib- 
rium. The oxygen-burning shell model is characterized 
by an upward kinetic energy flux, with a maximum am- 
plitude roughly equal to 5 % of the maximum enthalpy 
flux. Neutrino cooling is responsible for the shallow re- 



gion of negative luminosity. Nuclear burning is more im- 
portant and dominates Lnuc- A net heating results and 
the system is in thermal imbalance (Eq. 135] is not valid 
in this case). 

The bottom-left panels in Fig. [7] and Fig. [5] show 
the different terms of Eq. ([TO]) for models rg.3D.mr and 
ob.3D.mr. In the red giant case, we introduced in the 
r.h.s. the Newtonian cooling term q, see Eq. ([T]). This 
last term is only important at the top of the CZ, where it 
mimics the radiative cooling that would take place if we 
had a realistic photosphere. This convection is driven by 
cooling at the top. At the bottom of the convective zone, 
there is a rough balance between the radiative and the 
enthalpy flux. Radiation cools the gas in a shallow layer 
below the convective zone, and heats it at the bottom of 
the convective zone. The radiative cooling in a shallow 
layer below the convective zone is a different manifesta- 
tion of the bump seen in the radiative luminosity, this is 
discussed further in Sect. 14.61 The discussion of kinetic 
energy flux remains the same as in the previous section: 
kinetic energy is transported downward in the convective 
zone. The oxygen-burning shell model shows a rough bal- 
ance between nuclear heating and the divergence of the 
enthalpy flux: convection is driven by a heating from be- 
low. It should be noticed that the time derivative of the 
total energy contributes significantly to the balance, as a 
result from the thermal imbalance. Spurious oscillations 
are also present at the location of the steep gradients (at 
the convective boundaries), but we have checked that the 
consistency was good everywhere else. 

The bottom-right panels in Fig. [7] and Fig. |S]show the 
budget integrals over the convective zone for each model. 
In both case, a net heating term is present, in the red 
giant due to radiation, in the oxygen-burning shell due 
to nuclear burning. In the latter, this is counterbalanced 
by the evolution of the background state, through the 

change in time of the total energy and the — (p)Vr(Mr) 
term which is the rate of work by the mean pressure due 
to a global expansion/contraction of the star. In the red 
giant model, the radiative heating is mostly balanced by 
the cooling at the surface. The terms describing a global 
evolution of the background roughly balance each others, 
and describe how the background state is slowly evolving 
toward thermal equilibrium (see discussion in Sect. | 

4.2. The turbulent velocity field 
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Figure 10. Vortical correlation length-scales for the radial velocity, density and pressure fluctuations in model rg.3D.mr (left panel) and 
model ob.3D.mr (right panel). The dashed and dotted lines show the pressure and density scale-heights (respectively). 

In this section, we analyze the properties of the veloc- 
ity field, and consider relevant approximations for this 
field. We analyze in some detail the anisotropy of the 
Reynolds stresses. Our results motivate a two compo- 
nent decomposition to the flow, with plumes at large 
scales and isotropic turbulence at small scales. We show 
how the isotropic component is related to the observed ki- 
netic energy dissipation, consistent with the Kolmogorov 
cascade. 



i.e., the turbulent velocity field is solenoidal. Otherwise, 
for "deep" convection, where Lu ^ Hp, the appropriate 
approximation is: 



4.2.1. Approximations for deep and shallow convection 

Figures [3] and [H] show that the thermodynamical fluc- 
tuations in both stellar models are small relative to the 
background values. Therefore, for the purpose of analy- 
sis, we will often use simplified equations resulting from 
a linearization around the background state. In such a 
case, we shall use po, Pq, etc, instead of the notation 
(p), (P), etc, to emphasize that we are considering the 
linearized equations. 

The linearized continuity equation reads 



V • (pom') = 0. 
From this relation, we can deduce that 



(41) 



(42) 



The dilatation of the velocity field is due to the vertical 
motion of parcels in the b ackground stratification. 

iMeakin fc Arnet'^ ()20Q7[ ) compute the two-point corre- 
lation function of the radial velocity fluctuations: 



dtp' + V ■ {pm') = 0. 



{u',{t-r,6,4>)u'^{t;r + 6r,6,4>)) 

Urmsir)Urmsir + Sr) 



(43) 



(37) 



where 



— /„'2 



Both the red giant and oxygen-burning shell models 
are characterized by low Mach fiows. Low Mach turbu- 
lence is an inefficient producer of sound (Lighthill 1952). 
Therefore, we can neglect dtp' in the above equation and 
focus on fluctuatio ns haviiig a ti r ne scale longer than th e 
acoustic time scale (|Goughlll969l:lDutton fc Fichtlll 19690 . 
Furthermore, 



V • {pou') = poV/i • u',^ + poVru'r + Kdrpo, (38) 

where we have separated the horizontal flow (subscript 
h) from the vertical flow. The order of magnitude of the 
ratio of the third term to the second is: 



). The vertical correlation length- 



\u'^.drPo\ ^ Lu_ 
|poVr<| Hp' 



(39) 



where L„ is the characteristic vertical length-scale of ra- 
dial velocity perturbations, and Hp = —dr/d\npo is the 
density scale-height. For "shallow" convection, this ratio 
is small and the following approximation is justifled: 



scale is defined as the width at half maximum of the 
two-point correlation function . Figure [TU] shows 
the vertical correlation length-scale of the radial veloc- 
ity fluctuations (among others) computed at each radii 
in models rg.3D.mr and ob.3D.mr. The figure empha- 
sizes an important difference between the models: in 
the oxygen-burning shell model, the vertical correlation 
length-scale is everywhere smaller than the density scale- 
height, whereas in the red giant model it is larger in 
most of the convective zone. This suggests that the flow 
in the oxygen-burning shell model is much less affected 
by the stratification than the red giant model. There- 
fore, we will adopt the shallow convection approximation 
Eq. (Uni) to describe the turbulent velocity field in the 
oxygen-burning shell models. For the red giant models, 
the appropriate approximation is the one given by Eq. 

Equations (|40|) and (f4T|) are the basis for the Boussi- 
nesq and anelastic approximations. We stress that for 
both the oxygen-burning shell and the red giant models, 
these approximate models of the hydrodynamical equa- 
tions should not be used to model the fiow. For instance, 
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Figure 11. Eigenvalue analysis of the anisotropy tensor bijior model rg.3D.mr (left panel) and model ob.3d.mr (right panel). The figures 
show the radial profiles of the eigenvalues. The continuous line corresponds to the radial eigenvalue, the dashed lines correspond to the 
horizontal eigenvalues (see text). The two horizontal dashed lines indicate the range [—2/3, 2/3]. 



the computational domain in the oxygen-burning sheU 
models is not "shallow", i.e. (rout — rin)/Hp > 1, and 
in the red giant models the Mach numbers are too large 
near the surface. However, these approximations provide 
a useful theoretical framework for the analysis of the re- 
sults. 

4.2.2. Anisotropy of the Reynolds stresses 
We decompose the Reynolds stresses as 



(«) = -(fc)<5,, + (fc)6. 



(44) 



dominated by waves rather than turbulence. 

The anisotropy of the Reynolds stresses is another 
manifestation of the two components character of the 
flow, with convective plumes at large scales that dom- 
inate the vertical motion and an isotropic component 
at small scales. It is essential for the design of turbu- 
lence models to better understand the role of the different 
scales of t he flow, for instance in terms of the transpor t 
of energy (jCattaneo et al.lll991l : iBessolaz fc Brunll2011| ). 
We leave a detailed analysis for a future publication. 

Here, we outline an approach for decomposing the ve- 
locity field. We write: 



where (k) = is the specific turbulent kinetic 

energy and bij is a trace-less and symmetric tensor that 
characterizes the anisotropy of the Reynolds stresses. We 
extract the three eigenvalues (Ar,A0,A^) and eigenvec- 
tors of bij at each radii. It can be shown that, for any 
two eigenvalues a, (3, one has a > —2/3, /? > —2/3, and 
a-h/3<2/3. 

At each radius, one of the eigenvalue is unambigu- 
ously associated with a purely radial eigenvector, and the 
two other eigenvalues are roughly equal and associated 
with two purely horizontal eigenvectors. Therefore, the 
Reynolds stresses are axisymmetric, which is expected as 
the angular directions are homogeneous in the absence of 
rotation or magnetic field. The orientation of the hori- 
zontal eigenvectors is therefore not physically relevant. 
The radial profiles of the radial eigenvalue and of the 
horizontal eigenvalues Ag and A^ are shown in Fig. Illl for 
models rg.3D.mr and ob.3D.mr. These quantities describe 
the shape of the Reynolds stress tensor. In the bulk of 
the convective zone, where Ar > and Ag = A^ < 0, the 
shape is "rod-like" as the stress is maximum in the radial 
direction. Otherwise, when A^ < and Ae = A^ > 0, the 
shape is "disc-like". The radii at which the transition 
between the two shapes occurs (and where all eigenval- 
ues cancel) are ^ 2.2 x 10^^ cm and rout ^ 3.75 x 10^^ 
cm in the red giant model, and rjn — 0.45 x 10^ cm and 
''out ^ 0.85 x 10^ cm in the oxygen-burning shell model. 
At these radii, the vertical motions are defiected hori- 
zontally as they approach the convective boundaries. A 
second transition radius is seen at r = 2 x 10^^ cm in 
the red giant model and at r = 0.41 x 10^ cm in the 
oxygen-burning shell model, below this radius the flow is 



Ui = (ui) + v'l + w'l 



(45) 



where we have split the velocity fluctuation (■«■' in our 
usual notation) into a large scale, coherent structures 
component u", and an isotropic component w'l charac- 
terizing the small scales. We then assume that both com- 
ponents have zero averages and are not correlated: 



«)=0,Vz, 
(<<) =0'V(*,j) 



(46) 
(47) 
(48) 



With these assumptions, the velocity components cor- 
relation tensor can be split into two contributions: 



Furthermore, we have 



{w'/w'J) = ^hsoSij, 



(49) 

(50) 
(51) 



with (k) = /cpiumcs + fciso, where fcpiumcs and hso are 
the specific kinetic energy of the plumes and of the 
isotropic turbulence (respectively). The above hypoth- 
esis are not sufficient to determine uniquely the decom- 
position. Here, it w ill be e nough for our purpose to 
adopt the approach of lMeakin fc Arnett (,2007. ) : far from 
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the boundaries, we identify the horizontal flow with the 
isotropic component. Therefore, as a proxy, we can de- 
fine hso by 

2fciso-^ ((<) + (<)). (52) 
Furthermore, since (m^^) ~ {u'^), one has 

(wfw") W 2fcpiumcs'5il<5jl, (53) 

with 



2fcp,u,nc. = ^) ^"^'^^^"^^ (54) 

4.2.3. Kinetic energy damping 

The total kinetic energy dissipated per unit time is 
Ld ~ J {pf-d)dV, see Tables [T] and [2] The rate of dissipa- 
tion is not constant in time, but evolves with the flow. In 
the red giant model, we see a decrease in time of Ld to- 
ward the value quoted in the table as the model relaxes 
toward a quasi-steady state (see left panel of Fig. [2]). 
In the oxygen burning-shell model, Ld increases slowly 
with time as a result from the global heating due to the 
imbalan ce between nuclear h eating and neutrino cooling. 

As in lArnett et al.l (|2009( ) , we compute a dissipation 
length-scale Id and timescale as: 
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Figure 12. Comparison of the dissipation described by a ej. = 
u''^/ K law (dashed line) with the dissipation inferred from the data 
in model rg.3D.mr (continuous line). The dotted line shows the 
residual. 

the overall shape of the profiles agree very well. This 
suggests that a somewhat larger value of Id would give 
a better fit. The largest discrepancies are found close to 
the convective boundaries, where k^so is affected by the 
horizontal flow due to the interaction with the bound- 
aries. 

It should be emphasized that this approach works well 
in our case because we are in a statistically steady state. 
In situations where the kinetic energy changes rapidly, 
the dissipation rate lags behind the cascade rate given by 
u'^ /K due to the time needed to redistribute the k inetic 
energy over the inertial range (|Livescu et al.ll2009f ). 



j (pe^)dV^Mcz^ 



K,CZ 



1 Id 



(55) 

/ {ped)dV 2 Urms 

where Wrms is computed from the isotropic kinetic energy 
defined in the previous section, i.e. Urms = (2fciso)^^^- 
This yields Id = 7.7 x 10^^ cm in the red giant model, 
and Id ~ 0.4/cz ; and Id = 5.5 x 10^ cm in the oxygen 
burning-shell model, so Id ^ Icz- Furthermore, ^ 19 d 
in the red giant model and Td ^ 29 s in the oxygen- 
burning shell model. In both cases this is much shorter 
than the convectiv e turnover timescale. As discussed in 
lArnett et al.l (|2OO90 , this illustrates the strong dissipative 
character of turbulent convection. 

We now relate the kinetic energy dissipation inferred 
from the numerical simulations with global properties of 
the flow. The simplest approach is to use the formula 
for the rate of dissipation i n isotropic and homogeneous 
incompressible turbulence ()Popell2000l) : 



Ed 



,/3 



A 



(57) 



where u' and A are the rms v elocity and length-s cale 
of the most energetic eddies. lArnett et al.l (|2009f ) fit 
the dissipation in their oxygen-burning shell data with 
u' = (2A:iso)^^^ and by setting A to the dissipation length- 
scale Id derived above, see the right panel in their Fig. 2. 
In Fig. [121 we show that the same approach gives fairly 
good results for the red giant model as well. The am- 
plitude of the dissipation is slightly overestimated, but 



4.3. Magnitude of pressure fluctuations 

A comparison of Fig. [3] and Fig. |6] shows that p' /po ^ 
p' I po in the red giant model, whereas the oxygen-burning 
shell model is characterized by p'/po < p'/po- This is 
related to the background stratification. 

Within the framework of the approximations presented 
in Sect. 14.2.11 it is possible to obtain an elliptic equation 
for the pressure fluctuations (see Appendix IX)) : 

do' 

Ap' = -V : (pou' ® u') - g^. (58) 

The first term on the r.h.s. describes the generation of 
pressure fiuctuations by the Reynolds stresses, also called 
the "pseudo-sound" . The second term describes the gen- 
eration of pressure fluctuations by buoyancy effects. The 
relative orders of magnitude of the different terms can be 
written as: 



M 

Po 



Po 



Li 



HpLp Po 



(59) 



where Lp, L^, and Lp are typical length-scales for vari- 
ation of the pressure, velocity, and density fluctuations. 
We use the vertical correlation length-scales for u^, p' 
and p' shown in Fig. [10] to estimate the relative mag- 
nitude of these effects. We estimate that L-^/Ll^ ^ 0.2, 
Lp/ (HpLp) 2 in the red giant model, and L^/L'^ ^ 0.1, 
Lp/ {HpLp) ^ 0.4 in the oxygen-burning shell model. Al- 
though these estimates are only qualitative, they suggest 
that the pseudo-sound term is not the dominant source of 
pressure fluctuations in the convective zone (the dashed 
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Figure 13. Top panels: decomposition of the average mass flux (see text). Bottom panels: Splitting of the kinetic energy (thick lines) 
and acoustic (thin lines) fluxes into up-flow and down-flow components (the fluxes have been multiplied by 47rr^ and normalized by the 
luminosity). The continuous lines are the total fluxes. Left panels: model rg.3D.mr, right panels: model ob.3D.mr. 



line in Figs. [3] and [6] falls below the pressure curve after 
multiplication by Lp/L^). Therefore, pressure fluctua- 
tions are mainly due to buoyancy, and the order of mag- 
nitude analysis above shows that magnitude is related 
to the background stratification. If we loosely assume 
Lp ^ Lp (which is not quite correct in the red giant 
models), we have 



M 

Po 



Lp 

Hr, 



Po 



(60) 



which is in qualitative agreement with the numerical 
models. 

Although the order of magnitude estimates we carried 
out in this section have no predictive power, they suggest 
that Eq. ((58| provides a valuable framework to analyse 
pressure fiuctuations (iChassaing et al.ll2002[ ). We plan 
to push the analysis based on this equation further, with 
the aim of obtaining more quantitive predictions. 

4.4. The turbulent fluxes 
4.4.1. The turbulent mass flux 
The turbulent mass flux is defined as 



The top panels of Figure [13] illustrate this relation in 
models rg.3D.mr and ob.3D.mr. In the red giant model, 

one has (ur) ~ as the model is close to equilibrium. 

In this case, Eq. (|62|) implies that the mean fiow (ur) 
counter-balances the mass displaced by turbulence. The 
oxygen burning-shell model shows a quite significant ex- 
pansion of the background, i.e. (ur) > 0, driven by the 
imbalance between nuclear burning and neutrino cooling. 
There is still a mean fiow which tends to counter-act the 
effect of the turbulent mass flux. 

In the gravity field, the mass displaced by the turbu- 
lence induces a work which is one of the kinetic energy 
driving term introduced in Sect. 14.1.21 



Wb = {p'u' ■ g) = ~gf„ 



(63) 



/„, = {p'K). 



By definition, we have 



{p){Ur) = {pUr) = {p) {Ur) + (/<). 



(61) 



(62) 



where we took g outside of the averaging operator (the 
Cowling approximation). Note that energetically, Eq. 

(j62p implies that when {ur) ~ 0, the gravitational work 
done by the turbulence is canceled by the gravitational 
work done by the mean fiow, so that the total work done 
by gravity is zero. This is a direct consequence of mass 
conservation. 

4.4.2. Acoustic and kinetic energy fluxes 

Assuming that {p){Dt){ek) — 0, which is justified by 
the analysis in Sect. I4.1.2[ the integral version of the 
kinetic energy balance, Eq. (flS)) . reads 



16 



M. VlALLET ET AL. 




-0 5' 1 1 1 1 1 1 1 ' —0 015' 1 1 1 1 1 1 ' 

n.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 D.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 

i-(ciii) xlO'^ r(au) xlO' 

Figure 14. Comparison between the acoustic flux /p = (p'uj,) and pressure-dilatation Wp = (p'V ■ u') multiplied by the density scale- 
height Hp. Both terms were multiplied by Attt^ and normalized by the stellar luminosity. Left panel: model rg.3D.mr, right panel: model 
ob.3D.mr. 



^Trr^ifk + fp)= I {Wb + Wp~{ped))dV, (64) 

which simply says that transport by combined kinetic 
energy and acoustic fluxes is the residual between driv- 
ing and dissipation. It is not an explicit formula for 
the fluxes since the velocity field also enters in the r.h.s, 
both in the driving and in the dissipation. When /p and 
Wp are negligible , the e quation simplifies to Eq. (4) of 
iMeakin fc iWttl (|20Tof ). 

The bottom-left panel in Fig. [13] shows the radial pro- 
files of fk and fp in model rg.3D.mr. As discussed in Sect. 
I4T31 the kinetic energy flux fk is large and downward 
directed. The figure shows that fp is smaller, but not 
negligible, and upward directed. Furthermore, we show 
in the figure the splitting of fk and fp into the down- 
flow (uj, < 0) and upflow components [u'^ > 0). In both 
cases, the contribution from the downdrafts dominates 
significantly. This emphasizes the strong asymmetry of 
the flow which results from the large stratification. The 
figure shows that both components of the acoustic flux 
are positive, i.e. down-flows have p' < and up-flows 
have p' > 0. Since fp is not negligible and is upward di- 
rected (opposite to the kinetic energy flux) , a better un- 
derstanding of fp is necessary to understand what is set- 
ting the amplitude of the kinetic energy flux. The right 
panel in Fig. [T3] shows the kinetic energy and acoustic 
flux in model ob.3D.mr, which are both small compared 
to the enthalpy flux. The acoustic fiux shows a more 
complex behavior than in the red giant models. The ki- 
netic energy fiux is upward directed, and shows a large 
cancelation due to the approximate symmetry b etween 
upflows and downflows. IMeakin fc Arnet'B (|2Q1CI[ ) exper- 
imented with the oxygen burning shell models by chang- 
ing from a heating from below to a cooling from above, 
and found that the kinetic energy flux reversed direction 
(see their cl model). Qualitatively, this change of behav- 
ior is explained by the different direction of propagation 
of plumes, which propagate downward when triggered by 
cooling at the top, as in the red giant model. 

Finally, multiplying Eq. (l42t by p' and taking the av- 
erage one obtains: 



fp^{p'u';)=Hp{p'V ■u') = HpWp. (65) 

The left panel of Fig. [T3] shows that this relation holds 
to a very good degree in the red giant model, validat- 
ing a posteriori the use of the anelastic approximation 
(Sect. I4.2.ip . Equation ([65]) is very simflar to Eq. (l63l) : 
they both connect a flux, /,„ and /p, to a kinetic energy 
source term, Wb and Wp. For fp and Wp this is true 
only within the anelastic approximation. We will come 
back to this in Sect. 14.51 For the oxygen-burning shell 
model, consistently with the approximation V-u' = in- 
troduced in Sect. 14.2.11 we consider that Wp ~ which 
is a good approximation, fp is more complex than in the 
red giant model, where it is dominated by stratification 
effects. The right panel of Fig. [T4l shows that the effects 
of stratification, as described by the anelastic approxima- 
tion, reproduce the gross features of fp, but that other 
effects contribute. This could be due to pressure pertur- 
bations related to boundary effects, or to the background 
expansion, and requires further investigation. 

Note that we have the following relation: 



{?^V.{p,u'))=Wp~^, (66) 

which characterizes the deviation from the anelastic ap- 
proximation. 

4.4.3. Splitting of the enthalpy fiux 

The enthalpy flux describes the transport of heat by 
convection and is therefore important for stellar interior 
modeling. In this section, we use thermodynamical rela- 
tionships to study its relation to other turbulent fluxes. 
As two state variables are sufficient to determine the 
thermodynamic state, different expressions for the en- 
thalpy fiux can be derived. We choose the pressure as 
one of the variables, as pressure fiuctuations play a dif- 
ferent role in the two stellar models presented here. For 
the other variable, we will consider density, entropy, and 
temperature, respectively. For the oxygen-burning shell 
model, we describe composition effects in terms of the 
average number of nucleons A and free electrons Z per 
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Figure 15. Splitting of tlie enthalpy flux in models rg.3D.mr (left) and ob.3D.mr ( righ t). Top panels: decomposition with Eq. II68I I. 
Middle panels: decomposition with Eq. I I70I I. Bottom panels: decomposition with Eq. JTTJ. fh{Y') represents the sum of the composition 
terms due to A' and Z' . The difi'erent terms are multiplied by 47rr^, and normalized to the model luminosity. 



nucleu^. For the red giant model, these terms are zero 
because the composition was uniform. 
Formahy, the turbulent enthalpy flux is defined as 

fh = {p){h"u'^), but we have checked that it is iden- 
tical to fh = (p) {h'u[.) (the same holds for the other 
turbulent fluxes). Based on this second form, we first 
split the enthalpy flux in terms of density, pressure, and 
composition fluctuations. Assuming that fluctuations of 
these variables are small. 



® If the composition variables Y = 1/A and Yn = Z/A are used, 
most of the efl'ect is concentrated in the single variable Y. Ye is 
almost constant in the oxygen-burning shell. 



h' 



+ 



dh 
dp 
dh 
'dA 



p,A,Z 



A' 



dh 
dp 
dh 
' dZ 



P 



p,A,Z 



P,P,A 



(67) 



which leads to the following expression: 
P Fi 



fh — — 



P Ts 
de 
'dA 



1 



(P'<) + -^(P'<> 

1 a — -1- 



Oc 

P^P,z^^'<^'-PdZ 



P:P,A 



(Z'K), (68) 
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where the coefficients are evaluated using the background 
state. This relation illustrates how the enthalpy flux 
can be decomposed into separate contributions from the 
turbulent mass flux, the acoustic flux, and composition 
fluxes. The top panels in Figure [15] show how this de- 
composition compares with the numerical data. We find 
that decomposition holds to an excellent degree of 
approximation. In the red giant model, both the terms 
related to density and pressure fluctuations contributes, 
whereas in the oxygen-burning shell model, the term re- 
lated to density fluctuations provides the main contribu- 
tion. 

A second way to split the enthalpy flux is by introduc- 
ing the entropy flux in place of the turbulent mass flux. 
In this approach, we start from 



effects. This leads to a different expression of kinetic en- 
ergy driving, as shown below. 

In the linearized velocity equation, the acceleration on 
the r.h.s. is 



Po Po 
see Appendix El It can be also written as 

,P' , fp' Hpp' 



(72) 



(73) 



Po ^ Po Hp Pq 

(jBraginsky fc Roberts! [19951 ) . The physical meaning of 
the first term can be elucidated by computing its work: 



, , „ , 1 / dh 
h' =Ts' + -p' + — 
p oA 



s,p,z dZ 



(69) 



which, by the same arguments as above, leads to 



fh —TJs + fp 
dh 
OA 



P 



dh 
PdZ 



s.p.A 



{Z'K). (70) 



The middle panels in Fig. [15] shows how this relation 
compares with the numerical data. In both models, the 
agreement is very good. In the red giant model, the 
enthalpy flux mainly results from the entropy flux, with 
a non-negligible contribution from the acoustic flux. In 
the oxygen-burning shell models, the contributions from 
both the entropy and composition fluxes are important. 
Comparing this with the first splitting, it illustrates that 
density fluctuations are here both due to thermal effects 
(i.e. entropy fluctuations) and composition effects. 

Finally, it is useful to write the decomposition of the 
enthalpy flux by introducing temperature fluctuations, 
one obtains: 



-pou 



Po 



-V • {p'u') + • [pou' 



Po 



(74) 



This term gives rise to the transport by the acoustic 
flux and to a source term which is related to the deviation 
from the anelastic approximation (see Eq. [66]). This 
last term will contribute only when compressible effects 
become important, i.e. for Mg ^ 1. 

We identify the second term in Eq. ([73| with the buoy- 
ancy force, written in terms of density and pressure fluc- 
tuations. It can be written as: 



PO 



Ho Pi 







P_ 

PO 



FiPo 



P , 



PQCpQ dz 



where we used 



H„ 



TiHp 



6 dsQ 
Cp dz 



(75) 



(76) 



This term characterizes the deviation from adiabaticity 
of the background. We can make the connection with 
entropy and composition fluctuations, since we have the 
thermodynamical relationship 



h=pCp{T'u',) + {l-S)fp 



dh 



T,p,Z 



(A'K) 



dh 
"dZ 



T,p,A 



iZ'K). (71) 



where 5 — aT, with a the coefficient of thermal expan- 



sion at constant pressure. In the literature, pCp{T'u'^) 
is often considered as the enthalpy flux. Formally this 
is correct only for both (1) a polytropic gas [5 = 1) or 
pressure fluctuations are neglected [p' — 0), and (2) com- 
position effects are neglected. The left-bottom panel in 
Fig. [T51 shows that pCp{T'u'^) yields a very good approx- 
imation of Jh in the red giant model. The right-bottom 
panel illustrates that in the oxygen-burning shell model, 
composition effects are significant in this decomposition. 

4.5. Kinetic energy driving in compressible fluids and 
turbulent dissipation in the convective zone 

In Section l4.1.21 we expressed the driving of kinetic en- 
ergy in terms of Wb, related to density fluctuations, and 
Wp, related to pressure fluctuations. However, both den- 
sity and pressure fluctuations arise from different physi- 
cal processes: thermal effects (i.e. entropy fluctuations), 
dynamical effects (e.g. compressibility), and composition 
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Fi Po 



^-±s' 



s_ds_ 

Cp dyi 



Vi 



(77) 



with y[ = A, Z (with an implicit summation, yi — A, 
y2 ^ Z). 

Using these relations, the total kinetic energy driving 
Wi, A- Wp can be written as 



Wp + Wb =^Tf, 

Up 



Vad y ds 

Hp dyi 



P-.p 



6 djs) 
Cr, dz 



where we used 



Vac 



Thermal effects 
Composition effects 

Background effects 

Compressibility effects 
(78) 

The first two contributions 



in Eq. (1781) can be also written in terms of the turbu- 
lent mass flux and acoustic flux thanks to Eq. ([77]) . The 
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Figure 16. Splitting of the kinetic energy driving for model rg.3D.mr and ob.3d.mr. 



first three terms are related to buoyancy, which contrasts 
which the usual designation of Wb as the "buoyancy" 
driving. Figure \W\ illustrates this splitting of the driving 
for the numerical models. In our cases, both the devia- 
tion from adiabaticity and the compressibility effects are 
negligible. In the red giant model, there is no composi- 
tional effects and the above formulation is the most con- 
venient as it expresses the driving only in terms of the en- 
tropy flux. In the oxygen-burning shell model, composi- 
tion effects are important. However, since p'/Po < p' / Po, 
the driving can be expressed in terms of density fluctua- 
tions mainly. This is similar to the comment made with 
the right panels in Fig. [15] 

Finally, we can relate the enthalpy flux to the kinetic 
energy driving. In the oxygen-burning shell case, we have 
fromEq. ([55]): 



fh 



P Ti 



Hn 



(79) 



which connects kinetic energy driving (since Wp is negli- 
gible) with the enthalpy flux. For the red giant case, we 
have from Eq. ((78)) : 



V 



tip 

which we use Eq. ([70)) to obtain: 



(80) 



fh 



Vad 



When Wp is negligible, this gives Eq. ([70)) . 

In a quasi-steady state, we have Ld = j{Wb + Wp)dF, 
see Eq. ()32l) . Therefore, we can estimate Ld from 



Ld 



ad 



V 



fhdV 



dr 
cz Hp 



Vad-^cnw , 



(82) 



where Vad is the average value of the adiabatic gradi- 
ent, Lc is the average value of the enthalpy luminosity 
(47rr^//i) over the convective zone, and uh^, is the num- 
ber of pressure scale-heights in the convective zone. This 



overestimates Ld when Wp is not negligible, because of 
the factor Fa > 1 in Eq. (|8ip . Nevertheless, it shows 
that the turbulent dissipation is of the same order as 
the convective lumi nosity. This result was suggested by 
[Hewitt et al.l ([1975D . For instance, in the red giant mod- 
els, we have: Zc ~ 4 x 10^'' erg/s, Vad ~ 0.35, uh^ ~ 7.8, 
which gives ^ 11 x 10'^^ erg/s. In the oxygen-burning 
shell model, we have: Zc 5 x 10^^ erg/s, Vad 0.245, 
~ 2, which gives Ld ~ 2.45 x 10"^^ erg/s. Both val- 
ues agree well with the inferred values (see Tables [1] and 
[2]). As expected, the value obtained for the red giant is 
overestimated. 

4.6. Thermal effects and overshooting in the red giant 

model 

The Kclvin-Helmholtz timescale is defined as the ratio 
between the thermal energy and the luminosity of the 
star: 



tkh 



(83) 



L. ' 

where i^iut = / pe.idy is the total thermal energy. For the 
red giant the Kelvin- Helmholtz timescale is 2.1 x lO'^ yr. 
This is much longer than our simulations, which span 
roughly 8 years of model time. In fact we are not able to 
simulate over several thermal timescales, as it would be 
necessary to ensure that the models have reached ther- 
mal equilibrium. This is a common limitation to all nu- 
merical simulations of deep convective envelopes which 
include radiative cooling. One possible way to overcome 
this problem is to boost the luminosity, thereby bring- 
ing the dynam ical and thermal timescales closer to each 
other; e.g., D obler et al.l ([20061 ). A consequence is that 
the characteristic Mach number of the flow increases. 
The physical character of turbulent convection becomes 
very different, with properties closer to photospheric con- 
vection, as compressibility and superadiabatic effects be- 
come important. Furthermore, increasing the luminosity 
implies an increase in the thermal diffusivity, so that for 
a given Reynolds number, the Peclet number decreases. 
As discussed below, this will change the behavior at the 
convective boundaries. For these reasons, we prefer to 
use a realistic value of the luminosity. As shown in Sect. 
14. 1[ out of thermal equilibrium behavior can be taken 
into account in the framework of the mean-field equa- 
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Figure 17. Left panel: Dimensionless temperature gradients in model rg.3D.mr, the vertical dotted lines show the bottom convective 
boundary layer (see text). Right panel: entropy balance in model rg.3D.mr. 



tions. 

We characterize our red giant models with a global 
Peclet number, defined as 



Pe : 



Wrms^CZ 
X 



(84) 



where x is the average value of the radiative diffusivity 
over the convective zone, typically x = 9 x 10^^ cm^/s. 
We find Pe ~ 5200 in model rg.3D.mr. This large value 
characterizes a very efficient transport of heat, and thus 
a very efficient convection. 

In stellar evolution calculations, the locations of the 
convective zone boundaries are based on linear criteria 
for dynamical stability, namely the Schwarzschild or the 
Ledoux criteria (the latter takes composition gradients 
into account). How ever, due to inertia, fiuid parcels 
can cross this limit. iZahnI (|1991l ) prese nts an analyti- 
cal investigation of the problem (see also iSchmitt et al.l 
1 19841 : lRempelll200^ . Zahn describes as "penetrative con- 
vection" the process in which the superadiabatic region, 
grows in size due to an efficient thermodynamical mixing 
at the convective boundary. This is the case on the earth, 
where the planetary boundary layer grows in size during 
the day. Regarding t he connection of the unstable region 
to the stable interior, iZahnI ()1991t ) distinguishes between 
"overshooting" , in which the transition is made directly 
in a shallow thermal boundary layer, and "subadiabatic 
penetration" , in which the penetrative flow first estab- 
lishes a nearly adiaba t ic, ye t stable, region below the 
convective zone. IZahnI (|1991f) suggests that in stellar in- 
teriors, owing to large values of the Peclet number, the 
conditions for subadiabatic penetration are fulfilled. 

Pioneer ing numerical s t udies of the p roblem are pre- 
sented in lHurlburt_etIaD (il986L Il994f ). iHurlburt et"all 
()1994[ ) study how subadiabatic penetration/overshooting 
changes depending on the "stiffness" of the interface, 
a free parameter of their models. Their models show 
subadiabatic penetration for low values of the stiffness, 
whereas for large stiffness they only have an overshoot- 
ing layer. More recent 2D work bv Rogers fc Glatzmaier 
(2005) arrive to a similar conclusion iBrummell et al.1 
((200211 re visit the problem in 3D, using a similar ap- 
proach as IHurlburt et al.l ()1994[ ). None of their 3D mod- 
els show evidence for a subadiabatic region. The authors 
argue the reason is the lower filling factor of plumes in 



3D turbulent convection, resulting in lower local Peclet 
numbers. Surprisingly, their numerical models have a 
smaller convective region than in the initial model. 

How do the red giant models compare with these stud- 
ies ? We discuss here only the bottom boundary, as the 
analysis of the top boundary is undermined by our ar- 
tificial treatment of the surface cooling. The left panel 
of Fig. [T7] shows the dimensionless temperature gradi- 
ents in model rg.3D.mr. The stratification is very close 
to adiabatic in the bulk of the convective zone, owing 
to efficient convection. We show with two dotted lines 
the region where the enthalpy fiux has a negative bump 
(see left panel of Fig. jH]). We identify this region with 
the convective boundary layer. It has a radial extent of 
60% of the local pressure-scale height. The change in 
sign of the enthalpy fiux marks the start of the stably 
stratified (subadiabatic) region. The radius at which it 
happens is only slightly smaller than the location of the 
convective boundary in the initial model: our models do 
not show evidence for strong convective penetration. Do 
our models show evidence for subadiabatic penetration? 
Based on an inspection of the left panel in Fig. [IT] and 
of the profile of iV^, see right panel in Fig. [TJ we can 
see a shallow, nearly adiabatic region, which occupies 
roughly 30% of the convective boundary region. The 
thermal boundary layer, where the temperature gradi- 
ent connects smoothly to the radiative value, occupies 
th e remaining 70%. There fore, we obtain a similar result 
to IBrummell et al.l (|2002f ) : this model is characterized 
by overshooting. The size of our convective boundary 
layer is on the low side of the range of values result- 
ing from the parameter study of IBrummell et aTl (j2002i) . 
This suggests that our boundary is rather "stiff'. Note 
that this stiffness is a natural outcome of our models, 
and not an input parameter as in the above mentioned 
studies. Furthermore, we have a realistic thermal con- 
ductivity profile, depending on density and temperature 
rather than only on dep th. 

FoUowinglliS ()1991D . IBrummell et al.l (|2002D suggest 
that they would obtain subadiabatic penetration by mod- 
cling higher Peclet number fiows. The oxygen-burning 
shell models have a formally infinite Peclet numbeill[ 
as thermal diffusion is negligible, and give an insight 

^ Being negligible, radiative diffusion was removed from the code 
for increased efficiency. 
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Figure 18. Comparison of the kinetic energy balance for different resolutions. Top panels: models rg.3D.lr and rg.3D.mr. Bottom panels: 
models ob.3D.mr and ob.3D.hr. 



into the large Pcclet number limit. As discussed in de- 
tail in Meakin & Arnctt (200^, the oxygen-burning shell 
models show evidence for "turbulent entrainment" . The 
physical process is one in which the turbulent kinetic 
energy present at the boundary is converted to poten- 
tial energy as it draws material into the convection zone, 
mainly through shear instabilities and wave breaking. As 
a result, the stratification is weakened, and this leads a 
steady incr ease in the size o f the c onvective zone, see 
Fig. 4 in iMeakin fc ArnettI (|2007[ ). As such, this is 
the same effect as the convective penetration discussed 
above. However, convective penetration seems to be usu- 
ally thought as being due to the combined effect of large 
scale plumes, whereas turbulent entrainment describes 
the continuous erosion by the small scale turbulence at 
the interface. Both effects have the same signature: a 
negative buoyancy work {Wt < 0). In the first case, it 
results from buoyancy braking of the plumes; in the sec- 
ond case, it results from the work the flow is doing against 
gravity in the process of mixing the stable layer material. 
These effects are not easy to distinguis h in numerical sim- 
ulatio ns, and both may contribute. IMeakin &: Arnet^ 
(|2007f ) characterize the entrainment rate at convective 
boundaries based on the bulk Richardson number Ri^, 
which measures the stiffness of the interface by taking 
turbulence into account. This elucidates the behavior of 
the convective boundary at very large Peclet numbers. 

At lower Peclet numbers, how do non-adiabatic effects 
modify this process? The convective boundary layer in 
the red giant models is characterized by W^f, < and 
Wp < 0; see the middle-left panel of Fig. [T] The figure 



shows that Wb is dominant: kinetic energy is mostly con- 
verted into potential energy. Furthermore, the bottom- 
left panel of Fig. [7] shows that the divergence of the en- 
thalpy fiux leads to heating in the convective boundary 
layer. This is better analyzed in the framework of the en- 
tropy balance, Eq. (|17p. which is shown in the right panel 
of Fig. [T71 In the overshooting layer, the divergence 
of the entropy flux heats (note that in the overshoot- 
ing layer fh « T/^, see middle-left panel in Fig. [T5|) . 
However, it is compensated by cooling from radiation. 
As a consequence, a quasi-steady state in which non- 
adiabatic processes counter-balance the effects of turbu- 
lent entrainment is possible: the convective region does 
not increase in size as it does in the oxygen burning case 
where radiative effects are negligible. The bump in the 
radiative luminosity is another manifestation of this pro- 
cess. This can be seen also on the left panel in Fig. 
1171 where the convective boundary is characterized by 
a temperature gradient that is subadiabatic yet super- 
radiative: 



V,ad < V < V 



ad ■ 



(85) 



iZhang et al.l ()2012|) show how taking this effect into 
account improves the agreement of solar models w ith he- 
lioseismology data (see also Christensen-Dalsgaard et al.l 

Overshooting and subadiabatic penetration corre- 
spond to "thermally inhibited" turbulent entrain- 
ment/penetrative convection. Subadiabatic penetration 
requires large Peclet numbers. The estimate given by Eq. 
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Figure 19. Instantaneous snapshots of the radial velocity field in an horizontal plane in the middle of convective zone for models ob.3D.lr, 
ob.3D.mr, ob.3D.hr (from left to right). Light tones indicate upfiows, dark tones indicate downflows. Each subfigures shows the full angular 
domain 45° X 45°. 



(|M|) maybe misleading, as the scales which are involved 
in the subadiabatic penetration process can be much 
smaller, and characterized by lower "turbulent" Peclet 
numbers. Furthermore, the bulk Richardson number 
might be relevant for subadiabatic penetration as well. 
A better understanding of the relativ e effects of the spo- 
radic plumes (Mcakin &._Arnett 2007; IArnett et al.ll2009D 
that hit the convective boundary and of the continuous 
erosion by the turbulence is necessary. We suspect these 
concepts could shed ne w light on the results obtained in 
iBrummell et~all (|200l . 

4.7. Comparison of different resolutions 

We have performed the various analysis presented in 
the previous sections for different resolutions and found 
good agreement, as shown for instance by the various 
quantities summarized in Tables [T] and [2j We have not 
found any significant deviation in the physical results 
that could stem from resolution issues. The oxygen- 
burning shell models at the lowest resolution (ob.3D.lr) 
shows spurious oscillations in some averaged quantities. 
This seems to be related with difficulties Riemann based 
solvers have with stratification at low resolution. Even 
in this case, both the profiles and amplitudes of these 
quantities are actually in good agreement with the re- 
sults obtained from higher resolutions which are free of 
these problems. However, we do not have convergence 
in the mean fields in the narrow region of steep gradi- 
ents at the base of the oxygen-burning shell, which are 
unresolved in the models considered here. There, the 
dissipation and compositional mixing are affected at the 
grid scale by the numerical algorithm, which undermines 
the RANS analysis. The Riemann solver replaces a steep 
gradient by a contact discontinuity while the RANS anal- 
ysis requires continuity; this issue requires further study, 
although the general behavior is relatively sane and the 
discrepancy localized. 

Within the convective region, the mean-field analysis 
shows robust behavior regarding resolution. The most 
resolution sensitive diagnostic might be the kinetic en- 
ergy dissipation, which in our models is purely due to 
numerics at the grid scale. Figure [15] compares the ki- 
netic energy balance in models rg.3D.lr and rg.3D.mr, and 
models ob.3D.mr and ob.3D.hr. The balance looks nearly 
the same at different resolutions. The kinetic energy dis- 
sipation profiles are very similar, although the resolutions 
differ by a factor of two. This suggests that the kinetic 
energy dissipation is set by the large scale properties of 



the flow, and does not depends on the physics at small 
scales (here the grid scale). We interpret this as an in- 
dication that the dynamics in the turbulent convective 
zone is governed by the large scale dynamics, charac- 
terized by the coherent plumes which are well resolved 
even at our lowest resolution. This is consistent with 
the picture of the turbulent cascade (Richardson, If 9221 : 
iKolmogorovl I f 94 fi ) : at large Reynolds numbers the rate 
of dissipation is set by the energy injection at large scale, 
and is independent of the value of the viscosity. This is 
the so-called "dissipation anomaly". The viscosity sets 
the scale at which dissipation occurs, here the grid scale. 

The turbulent regime obtained in numerical simula- 
tions (ours and others) is characterized by coherent 
plumes which propagate (upward or downward) over a 
significant fraction of the convective region, if not the 
whole region. They govern the large scale dynamics and 
are key in guiding the modeling of the highly non-local 
and non-isotropic transport properties of the fiow (see 
e.g. Rcmpcl 2004; LcsafFre et al. 20 05t iBelkacem et al.l 
i2006:,Kupka fc Robinson,.2007£ .Meakin fc Arnettll2007D . 
However, one should bear in mind that our numerical 
simulations have non-dimensional numbers (e.g.. Re, Pr, 
Ra) which are orders of magnitude different from the val- 
ues relevant to stellar hydrodynamics. Figure 1191 shows 
snapshots of the flow in the oxygen-burning shell model 
for three different resolutions. Although the flow is char- 
acterized by structures at smaller and smaller scales, we 
do not see any evidence for a different global behavior of 
the flow. Whether a transition to a different regime oc- 
curs at (much) larger resolution is an outstanding prob- 
lem. A similar question arises in the study of the simpler, 
but not less fundamental, Rayleigh-Benard convection 
problem. The quest for an understanding and charac- 
terization of the "ultimate" state of turbulent Rayleigh- 
Benard convection is the focus of muchjjxpcrimcntal and 
theoretical work (see reviews bv [Siggiail994i:IAhlers et all 
l2009f ) . Although there are significant physical differences 
with the stellar case, e.g., stemming from the different 
nature of boundaries, the extremely low values of the 
Prandtl number, or the effects of compressibility, it can 
be expected that a better understanding of turbulent 
Rayleigh-Benard convection will provide valuable insight 
into the turbulent regime at which stellar c onvection op- 
erates (see e.g. discussion in lSpruitlll997[ ). Theoretical 
studies support t he existence of these plumes in stellar 
convective zones (jSimon fc Weisslll99lHRieutord fc ZahnI 
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|1995[) . Numerical models of the propagation of a plume 
through the adiabatic stratifica tion show the develop- 
ment of secondary instabilities (jRasti 119981 : iClvne et al.l 
120071 ). but they deal with an idealized situation as the 
interaction between plumes clearly dominates in our nu- 
merical models. Observationally, these plumes are too 
deep and too small to be de tected by current heli oseis- 
mology measurements; see (jHanasoge et al.l I2012D and 
references therein. 

5. CONCLUSION 

This paper presented 3D models of the turbulent con- 
vection in the envelope of a red giant star and in the 
oxygen-burning shell of a supernova progenitor. The two 
models differ significantly in their physical properties: 
they have radically different equations of state, the effects 
of thermal diffusion is negligible in one but important in 
the other, one model is multi-fluid and includes nuclear 
burning, whereas the other is mono-fluid and has cooling 
at the surface. Their common point, which is the focus 
of this work, is the presence of a turbulent convective 
zone which dynamics is controlled by the hydrodynami- 
cal equations. Finally, two different numerical methods 
and codes were used to produce the data. To deal with 
such a heterogenous set of data, we developed in Sect. |3] 
a set of ID horizontally-averaged equations that provides 
a framework for a systematic analysis of hydrodynamical 
simulations. We showed in Sect. 14.11 that our numerical 
models show good consistency with the physical equa- 
tions, although we identified spurious effects localized in 
the region of the steep, unresolved gradients present in 
the oxygen-burning model. 

Both our models are characterized by low Mach flows, 
so that compressible effects are negligible. Similarly, 
both models have large Peclet number (formally infinite 
in the oxygen-burning shell case) so that the deviation 
from adiabaticity is small in the convective zone, and has 
no effect on the dynamicfl However, our analysis is gen- 
eral and applies also to fiows with larger Mach number 
and superadiabatic stratifications, two conditions which 
are found in photospheric convective regions. We plan to 
apply a similar analysis to photospheric convection. Our 
mean-field analysis emphasized very similar behavior in 
both stellar models, without noticeable dependance on 
the numerical resolution. Both the radial and kinetic en- 
ergy balances are in a statistically steady state, whereas 
the total energy balance shows an evolution of the back- 
ground on a longer timescale. This is due to the clear 
separation between the dynamical timescale of the sys- 
tem (the turnover timescale) and the nuclear/thermal 
timescale of the models. The kinetic energy dynamics 
can be understood as a balance between driving at large 
scales and dissipation at small scales, connected by the 
turbulent cascade. Understanding the spatial distribu- 
tion of the driving and of the dissipation is important for 
insight into the non-locality of convection. 

We have shown that the differences between the red gi- 
ant and oxygen-burning shell models stem mainly from 
the degree of stratification. This led us to introduce the 
distinction between "shallow" convection, in which the 
velocity correlation length-scales are less than the den- 

* In the red giant model, this would not be the case near the 
surface if it was modeled realistically. 




Figure 20. Global energy balance in turbulent convection. The 
arrow directions are consistent with the labels, i.e. for positive val- 
ues the energy flows along the arrows. Continuous lines emphasize 
the reversible character of the processes, whereas the dashed line 
denotes the irreversible character of kinetic energy dissipation. 

sity scale-height, and "deep" convection in which they 
are of the same order, or larger. As shown in Sect. 14.31 
the effect of stratification on the dynamics can be under- 
stood in terms of the magnitude of pressure fluctuations. 
We discussed how this impacts the mean-field balances: 
for deep convection pressure-dilatation becomes a non- 
negligible source of kinetic energy, and the acoustic flux 
contributes to the transport of kinetic energy and en- 
thalpy. We showed in Sect. 14.51 the connection between 
the transport of enthalpy, the rate of production of tur- 
bulent kinetic energy, and flnally the rate of turbulent 
dissipation at small scales. This should not be surprising: 
large scale transport of enthalpy needs motion, motion 
becomes turbulent, and turbulence dissipates kinetic en- 
ergy at small scales. As a consequence, we flnd that in 
a quasi-steady state the rate of dissipation of turbulent 
kinetic energy is of the order of the convective luminos- 
ity. This is consistent with the rate of damping of ki- 
netic energy inferred from our models. How turbulent 
dissipation affects stellar evolution is an open question. 
Figure [20] puts the role of turbulent dissipation in per- 
spective regarding the global conservation of energy. The 
convective instability taps energy of the unstable stratifi- 
cation and converts potential energy into kinetic energy, 
the work done by the background flow {ur) on the mean 
background (p) converts internal energy into potential 
energy, and turbulent dissipation closes the loop by con- 
verting kinetic energy into internal energy. In a statis- 
tically steady state, the amount energy per unit time 
which is "flowing" in these channels is Wb (see Eqs. |32] 
and [621) ■ 

We have discussed in Sect. 14.61 the overshooting pro- 
cess which is observed in the red giant models. Compar- 
ing with the oxygen-burning models, we discussed how 
turbulent entrainment relates to the more classical con- 
cepts of overshooting and subadiabatic penetration used 
in the stellar context. Further investigations on the rel- 
ative effect of the plumes and turbulent entrainment is 
desirable. 

A turbulent model for stellar convection should be able 
to reproduce the behavior of the mean-field equations 
without the need to resort to expensive 3D simulations. 
In the RANS framework, evolution equations can be de- 
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rived for the turbulent fluxes. The resulting equations 
involve higher order terms, for which additional evolu- 
tion equations can be derived. This leads to a hierarchy 
of equations which have to be closed with appropriate 
relations. Future work will apply the same systematic 
study to higher order equations, aiming at identifying 
the most important terms to guide the closure strategy 
(Mocak et al., in preparation). We intend to release pub- 
licljl3 our RANS analyzed data and provide analysis and 
plotting subroutines as open source materials. 

We plan to extend this work by including rotation and 
magnetic field. As mentioned in Sect. [3l the formula- 
tion of suitable ID mean-field equations is not possible 
when rotation or magnetic field are included. Neverthe- 
less, if axisymmetric instabilities are not important, it 
is possible to average the equations over the azimuthal 
direction, resulting in a set of 2D mean-field equations. 
This opens the possibility of performing stellar evolution 
in 2D, with an appropriate treatment of these effects (in 
the lim it of low r otation rate due to sp herical geometry) , 
see e.g. iDeupred (|I990, 20011: lLi et al . (2006, 2009). The 
MUSIC code, which is based on time-implicit methods, 
provides the ideal framework for that. 
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appendix 

elliptic equation for the pressure 

We start from the momentum equation: 



Po ( dtu + u • Vu ) — — Vp' 



p'g, 



(Al) 



where we have neglected density fluctuations in front of the Lagrangian derivative, and we have subtracted the 
hydrostatic background. Taking the divergence of this equation removes the time derivative both in the Boussinesq 
(V • u=Qi) and in the anelastic (V • {pou)=Gi) approximations: 



where we considered that g was constant. We write this equation as 



or 



dp' 



(A2) 



(A3) 
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